1. Загрузите датасет life_expectancy_data.RDS (лежит в папке домашнего задания). Это данные с основными показателями, через которые высчитывается ожидаемая продолжительности жизни по метрике World Development Indicator на уровне стран2. В данных оставлены строки, относящиеся к положению женщин в 2019 г.

df <- readRDS('life_expectancy_data.RDS')
head(df)
##                Country Year Gender Life expectancy Unemployment
## 1:         Afghanistan 2019 Female          66.388    14.065000
## 2:             Albania 2019 Female          80.201    11.322000
## 3:             Algeria 2019 Female          78.133    18.629999
## 4:              Angola 2019 Female          64.039     7.835000
## 5: Antigua and Barbuda 2019 Female          78.104     8.256285
## 6:           Argentina 2019 Female          79.995    10.698000
##    Infant Mortality          GDP          GNI
## 1:             42.9  18799450743  19098305631
## 2:              7.7  15400242875  15198661275
## 3:             18.6 171767403748 167568133680
## 4:             44.5  89417190341  81900890781
## 5:              5.1   1687533333   1581477852
## 6:              7.0 451932356086 434037404778
##    Clean fuels and cooking technologies Per Capita
## 1:                                 36.0   494.1793
## 2:                                 80.7  5395.6595
## 3:                                 99.3  3989.6683
## 4:                                 49.6  2809.6261
## 5:                                100.0 17376.6497
## 6:                                 99.8 10056.6379
##    Mortality caused by road traffic injury Tuberculosis Incidence
## 1:                                    15.9                    189
## 2:                                    11.7                     16
## 3:                                    20.9                     61
## 4:                                    26.1                    351
## 5:                                     0.0                      0
## 6:                                    14.1                     29
##    DPT Immunization HepB3 Immunization Measles Immunization Hospital beds
## 1:               66                 66                   64     0.4322222
## 2:               99                 99                   95     3.0523077
## 3:               91                 91                   80     1.8000000
## 4:               57                 53                   51     0.8000000
## 5:               95                 99                   93     2.5806250
## 6:               86                 86                   94     4.6100000
##    Basic sanitation services Tuberculosis treatment Urban population
## 1:                  49.00617               91.00000           25.754
## 2:                  99.18307               88.00000           61.229
## 3:                  86.13850               86.00000           73.189
## 4:                  51.39362               69.00000           66.177
## 5:                  85.52290               72.26667           24.506
## 6:                  91.40773               47.00000           91.991
##    Rural population Non-communicable Mortality Sucide Rate continent
## 1:           74.246                       36.2         3.6      Asia
## 2:           38.771                        6.0         2.7    Europe
## 3:           26.811                       12.8         1.8    Africa
## 4:           33.823                       19.4         2.3    Africa
## 5:           75.494                       17.6         0.8  Americas
## 6:            8.009                       12.1         3.3  Americas
colnames(df) <- gsub(" ", "_", colnames(df))
summary(df)
##    Country               Year         Gender          Life_expectancy
##  Length:195         Min.   :2019   Length:195         Min.   :55.49  
##  Class :character   1st Qu.:2019   Class :character   1st Qu.:70.02  
##  Mode  :character   Median :2019   Mode  :character   Median :77.55  
##                     Mean   :2019                      Mean   :75.52  
##                     3rd Qu.:2019                      3rd Qu.:80.95  
##                     Max.   :2019                      Max.   :88.10  
##   Unemployment    Infant_Mortality      GDP                 GNI           
##  Min.   : 0.178   Min.   : 1.40    Min.   :1.884e+08   Min.   :3.754e+08  
##  1st Qu.: 3.735   1st Qu.: 5.35    1st Qu.:1.117e+10   1st Qu.:1.094e+10  
##  Median : 5.960   Median :13.50    Median :3.967e+10   Median :4.009e+10  
##  Mean   : 8.597   Mean   :19.61    Mean   :4.660e+11   Mean   :4.864e+11  
##  3rd Qu.:10.958   3rd Qu.:30.23    3rd Qu.:2.476e+11   3rd Qu.:2.457e+11  
##  Max.   :36.442   Max.   :75.80    Max.   :2.143e+13   Max.   :2.171e+13  
##  Clean_fuels_and_cooking_technologies   Per_Capita      
##  Min.   :  0.00                       Min.   :   228.2  
##  1st Qu.: 34.50                       1st Qu.:  2165.3  
##  Median : 80.70                       Median :  6624.8  
##  Mean   : 65.98                       Mean   : 16821.0  
##  3rd Qu.:100.00                       3rd Qu.: 19439.7  
##  Max.   :100.00                       Max.   :175813.9  
##  Mortality_caused_by_road_traffic_injury Tuberculosis_Incidence
##  Min.   : 0.00                           Min.   :  0.0         
##  1st Qu.: 8.20                           1st Qu.: 12.0         
##  Median :16.00                           Median : 46.0         
##  Mean   :17.06                           Mean   :103.8         
##  3rd Qu.:24.00                           3rd Qu.:138.5         
##  Max.   :64.60                           Max.   :654.0         
##  DPT_Immunization HepB3_Immunization Measles_Immunization Hospital_beds   
##  Min.   :35.00    Min.   :35.00      Min.   :37.00        Min.   : 0.200  
##  1st Qu.:85.69    1st Qu.:81.31      1st Qu.:84.85        1st Qu.: 1.301  
##  Median :92.00    Median :91.00      Median :92.00        Median : 2.570  
##  Mean   :87.99    Mean   :86.76      Mean   :87.31        Mean   : 2.997  
##  3rd Qu.:97.00    3rd Qu.:96.00      3rd Qu.:96.50        3rd Qu.: 3.773  
##  Max.   :99.00    Max.   :99.00      Max.   :99.00        Max.   :13.710  
##  Basic_sanitation_services Tuberculosis_treatment Urban_population
##  Min.   :  8.632           Min.   :  0.00         Min.   : 13.25  
##  1st Qu.: 62.919           1st Qu.: 73.00         1st Qu.: 41.92  
##  Median : 91.144           Median : 82.00         Median : 58.76  
##  Mean   : 77.380           Mean   : 77.57         Mean   : 59.12  
##  3rd Qu.: 98.582           3rd Qu.: 88.00         3rd Qu.: 78.02  
##  Max.   :100.000           Max.   :100.00         Max.   :100.00  
##  Rural_population Non-communicable_Mortality  Sucide_Rate        continent 
##  Min.   : 0.00    Min.   : 4.40              Min.   : 0.300   Africa  :52  
##  1st Qu.:21.98    1st Qu.:11.85              1st Qu.: 2.050   Americas:38  
##  Median :41.24    Median :17.20              Median : 3.500   Asia    :42  
##  Mean   :40.88    Mean   :17.05              Mean   : 4.802   Europe  :48  
##  3rd Qu.:58.08    3rd Qu.:22.10              3rd Qu.: 6.600   Oceania :15  
##  Max.   :86.75    Max.   :43.70              Max.   :30.100

2. Сделайте интерактивный plotly график любых двух нумерических колонок. Раскрасть по колонке континента, на котором расположена страна

fig <- plot_ly(data = df, x = ~Life_expectancy, y = ~Urban_population,color=~continent, colors = "Set1")
fig
## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter
## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

3. Проведите тест, на сравнение распределений колонки Life expectancy между группами стран Африки и Америки. Вид статистического теста определите самостоятельно. Визуализируйте результат через библиотеку rstatix.

ggplot(df, aes(x=Life_expectancy, color = 'Life_expectancy'))+  
  geom_histogram(aes(y = ..density..), bins = 7,  fill = '#67B7D1', alpha = 0.5) +  
  geom_density(color = '#67B7D1') +  
  geom_rug(color = '#67B7D1') + theme_classic() +
  scale_color_manual(values = c('density' = '#67B7D1'))+labs(title="Distplot for Life expectancy",
        x ="LE (years)", y = "Density")

library(rstatix)
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
## 
##     filter
qqnorm(df$Life_expectancy, pch = 1, frame = FALSE)
qqline(df$Life_expectancy, col = "steelblue", lwd = 2)

print(df%>%shapiro_test(Life_expectancy))
## # A tibble: 1 × 3
##   variable        statistic           p
##   <chr>               <dbl>       <dbl>
## 1 Life_expectancy     0.944 0.000000712
# Несмотря на то, что отклонения от нормальности незначительные визуально (кроме обрезанного распределения, что и на qqплоте съезжает), тест на нормальность переменная не проходит. Т.к. это дз, я возьму MWU тест, но в реальной жизни, мне кажется, t test можно было бы тут использовать
library(ggpubr)


stat.test <- df %>% filter(continent %in% c("Africa", "Americas")) %>% 
  wilcox_test(Life_expectancy ~ continent, paired = FALSE)

p <- ggboxplot(
  df %>% filter(continent %in% c("Africa", "Americas")), x = "continent", y = "Life_expectancy", 
  color = "continent", palette =c("#00AFBB", "#FC4E07"),  add = "jitter"
)
p +stat_pvalue_manual(stat.test, label = "T-test, p = {p}", 
                      y.position = 85)

# Не совсем поняла, что нужно было сделать, надеюсь, сделала верно. Отличие стат значимое, как видим, на уровне 0.05

4. Сделайте новый датафрейм, в котором оставите все численные колонки кроме Year. Сделайте корреляционный анализ этих данных. Постройте два любых типа графиков для визуализации корреляций.

library(corrplot)
## corrplot 0.92 loaded
df_numeric <- df %>% 
  select_if(is.numeric) %>% 
  select(-Year)

correlation_matrix <- cor(df_numeric, use = "complete.obs")
# Баблплот
corrplot(correlation_matrix, method="circle", type="lower", 
         tl.col="black", tl.srt=45, diag=FALSE, 
         col = COL2('PiYG'), tl.cex = 0.4)

# Хитмап
corrplot(correlation_matrix, method="color", 
         tl.col="black", tl.srt=45, diag=FALSE, 
         col = COL2('PiYG'), tl.cex = 0.4)

### 5. Постройте иерархическую кластеризацию на этом датафрейме.

library(stats)
library(gplots)
## 
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
## 
##     lowess
set.seed(42)
df_scaled <- scale(df_numeric)
d <- dist(df_scaled)
hc <- hclust(d, method = "complete")
plot(hc, cex = 0.6, hang = -1)

### 6. Сделайте одновременный график heatmap и иерархической кластеризации. Содержательно интерпретируйте результат.

mat_clipped <- pmax(pmin(as.matrix(df_scaled), 4), -4)
row.names(mat_clipped) = df$Country

heatmap.2(mat_clipped, Rowv = as.dendrogram(hc), 
          symm = FALSE, trace = "none", 
          col = bluered(20), 
          margin = c(10, 10), 
          main = "Heatmap with Hierarchical Clustering", cexRow=0.3)

#### Интерпретация

Признаки кластеризуются по 2-м группам: 1) Признаки, связанные со смертностью или в дальнейшем “смертонстные” - все виды смертности попали сюда. Также сюда попали, что отчасти логично, безработица, заболеваемость туберкулезом, лечение туберкулеза и популяция деревни. Эти признаки “горят” в странах Африки и почти все понижены в ряде стран Европы. Интересно, что в ряде стран Европы сильно понижено лечение туберкулеза: это, возможно, говорит об очень низкой заболеваемости в том числе. 2) Признаки “комфорта” в свою очередь тоже бьются на 3 группы: - Признаки “иммунизации”, которые почти у всех стран более-мене на одном уровне, кроме ряда стран Африки, Азии в самом низу плота, у которых они явно понижены - “Городские” признаки: городская популяция, технологии и санитарные сервисы, а также, что интересно, продолжительности жизни. Видимо, эти признаки говорят о доступности базовых удобств, таких как чистая вода, и потому попали с продолжительностью жизни в один кластер. Эти фичи понижены в странах Африки и ряде стран Азии и повышены у всех остальных. - “Экономические”: ВВП, ВНД и ВВП на душу населения, а также риск суицида и кол-во мест в больницах – достаточно контринтуитивно попали сюда. Здесь почти все страны равны, помимо ряда, видимо, крайне развитых экономик Европы, Азии, Америки, в которых признаки повышены. Риск суицида достаточно шумный признак, что говорит, возможно, о том, что причины для суицида могут быть в любой экономике. А вот кол-во койкомест в больницах повышен в ряде стран Европы, при этом не в тех же, в которых явно самые высокие экономические скоры. Возможно, это постсовесткие страны. Страны в свою очередь группируются на 6: Группа 1: “смертностная” группа признаков повышена, “иммунизация” на нормальном уровне, “городские” и “экономические” признаки снижены или на среднем уровне, в группу попало множество стран Африки. Группа 2: повышена популяция деревни из “смертностных” признаков, остальные на среднем уровне, повышены “городские” признаки, “экономические” признаки снижены или на среднем уровне. Сюда попали страны Азии, Европы, Америк, Океании Группа 3: понижены “смертностные” признаки, повышены “городские” и на среднем уровне “экономические”. Здесь бол-во стран, видимо, из Америки Группа 4: самые низкие “смертностные”, высокие “городские” и самые высокие “экономические” - в основном страны Европы, а также немного из Америки и Азии. Группа 5: “смертностные” признаки самые высокие, самая низкая “иммунизация”, низкие “городские” признаки и средние “экономические”. Сюда попали страны Африки, такие как Ангола, Конго, Южный Судан. При желании эта группа вполне себе вписывается в первую, единственное отличие это “иммунизация” и более плохие показатели в целом. Группа 6: 2 аутлаера и один из них Китай. Одни из самых высоких “экономиечских” характеристик, из-за ВВП и ВНД, но не ВВП на душу населения. Остальные признаки на среднем уровне.

7. Проведите PCA анализ на этих данных. Проинтерпретируйте результат.

row.names(df_scaled) = df$Country
pca_result <- prcomp(df_scaled, center = TRUE, scale. = TRUE)
summary(pca_result)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6    PC7
## Standard deviation     2.7526 1.4841 1.3952 1.17177 1.08375 0.96347 0.9288
## Proportion of Variance 0.3988 0.1159 0.1025 0.07227 0.06182 0.04886 0.0454
## Cumulative Proportion  0.3988 0.5147 0.6172 0.68945 0.75126 0.80012 0.8455
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.85740 0.69263 0.68937 0.59106 0.54986 0.47085 0.36596
## Proportion of Variance 0.03869 0.02525 0.02501 0.01839 0.01591 0.01167 0.00705
## Cumulative Proportion  0.88421 0.90946 0.93447 0.95286 0.96877 0.98044 0.98749
##                           PC15    PC16    PC17    PC18      PC19
## Standard deviation     0.34546 0.26941 0.20224 0.06968 2.876e-16
## Proportion of Variance 0.00628 0.00382 0.00215 0.00026 0.000e+00
## Cumulative Proportion  0.99377 0.99759 0.99974 1.00000 1.000e+00

Интерпретация

Видим, что данные достаточно шумные, в первой компоненте всего 40% объясненной дисперсии, хотя 75% объясненной дисперсии накапливаются к всего лишь 5-й компоненте.

plot(pca_result, type = "lines")

По графику локтя перелом происходит уже на 2-й компоненте.

library(factoextra)
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
fviz_pca_biplot(pca_result,
                label="var",
                habillage = df$continent)+theme_classic()

#### Интерпретация

Видим, что первые две компоненты PCA, хотя не очень сильно информативные (в первой всего 40% объясненной дисперсии), похожим образом группируют характеристики: иммунизационные признаки формируют вектора с одним направлением, как и смертности, кроме суицида. Страны Африки, в основном группируются в направдении повышения “смертностных” признаков и в сторону от “иммунизации”. Страны Европы группируются в месте роста “экономических”. Кстати, что интересно, PCA более логично сгруппировал количество больничных мест вместе с другими признаками “комфорта” - наличием очистительных технологий и санитарии. Америки группируются в центре, т.е. не совсем объясняются компонентами. Здесь “аутлаером” является США, точка в крайне правой верхней части графика. Одна из стран Америки, также отстоит от других в левой верхней части. Страны Азии разнесены по первой компоненте и не образуют определенную группу. То же можно сказать об Океании, она находится посередине с тильтом влево вниз (по вектору деревенского населения). Но эти гипотезы мы также в дальнейшем подтвердим или опровергнем, построив биплот в плотли.

8. Постройте biplot график для PCA. Раскрасьте его по значениям континентов. Переведите его в plotly. Желательно, чтобы при наведении на точку, вы могли видеть название страны.

loadings <- pca_result$rotation[, 1:2]
loadings <- data.frame(Variable = rownames(loadings), loadings)
scaling_factor <- 3  #для лучшей визуализации

biplot_data <- data.frame(PC1 = pca_result$x[,1], 
                          PC2 = pca_result$x[,2], 
                          Continent = df$continent,
                          Country = df$Country)

p <- ggplot() +
  geom_point(data = biplot_data, aes(x = PC1, y = PC2, color = Continent, text = Country)) +
  geom_segment(data = loadings, aes(x = 0, y = 0, xend = PC1 * scaling_factor, yend = PC2 * scaling_factor), arrow = arrow(type = "closed", length = unit(0.2, "inches")), color = "black") +
  geom_text(data = loadings, aes(x = PC1 * scaling_factor, y = PC2 * scaling_factor, label = Variable), hjust = 1.1, vjust = 1.1) +
  theme_classic()

ggplotly(p, tooltip = "text")

9. Дайте содержательную интерпретацию PCA анализу.

Интерпретация

Вооружившись оптикой плотли, проанализируем результат с учетом выдвинутых ранее гипотез: 1. Страны Африки - Страны из самого “плохого” кластера преимущественно концентрируются в верхнем левом углу - Североафриканские страны в основном концентрируются в центре - Другие образуют облако слева внизу, не так далеко слева, как первые 2. Cтраны Америки - 2 аутлаера: США, самая правая сверху, и Гаити, слева сверху. - Остальные, в основном, концентрируются всередине, с небольшим сдвигом вправо. Интересно, что одни из самых бедных стран Америки, Боливия и Венесуэла, сдвинуты влево вверх, но не отстают слишком далеко от других, как например Гаити. 3. Страны Азии - В основном образуют единое облако, разнесенное по оси X: самая правая точка - Япония, левая - Афганистан. Китай не является явным аутлаером, как по результатам кластеризации. В целом можно разделить страны на: - более “центровые” (н-р Малазия, Ирак, Иордан, КНДР, Индия), - больше “справа сверху” (Япония, Китай, Корея, Сингапур – практически группа дальневосточной Азии, кроме КНДР), - больше “слева сверху” вместе с группой 5 Африки (Сирия, Филиппины, Йемен, Пакистан, Афганистан – крайне разнородная группа, но отчасти объединена наличием недавних войн, кроме Филиппин и Пакистана - хотя последняя известна потухшим индо-пакистанским конфликтом) - и больше “слева снизу” (Мьянма, Непал, Камбоджа, Бангладеш - буддисткие южноазиатские страны, крайне бедные). 4. Страны Европы - Выделяется группа концентрирумая справа сверху - страны западной Европы, а также часть пост-советских, например Россия, страны Балтии, Беларусь - Также облако чуть правее центра: страны пост-советские или из советского блока: Румыния, Украина, сербия, Грузия, Армения, Северная Македония. Видимо, чуть более бедные страны с теми же тенденциями, что в первой группе, т.к. как будто они идут немного параллельно. - “Маленькие” страны: неожиданно, Лихтенштейн, Косово – регион располагающийся наиболее слева среди европейских, Черногория, острова: Фареры и Нормандские, а также Босния и Герцеговина. 5. Страны Океании - Аутлаер Папуа Новая Гвинея - Остальные деляться на группы: - наиболее справа - Австралия, Новая Зеландия, и неожиданно Бруней, хотя предыдущие значительно ближе друг к другу; - наиболее сверху: Гуам, Новая Каледония, Французская полинезия - территории Франции или США в океании; - остальные находятся слева и разнесены по оси Y

10. Сравните результаты отображения точек между алгоритмами PCA и UMAP.

library(umap)

umap_result <- umap(df_scaled)
df_umap <- data.frame(UMAP1 = umap_result$layout[,1], UMAP2 = umap_result$layout[,2], continent = df$continent, Country=df$Country)
p_umap <- ggplot(df_umap, aes(x = UMAP1, y = UMAP2, color = continent, text = Country)) +
geom_point() +
ggtitle("UMAP")+theme_classic()


ggplotly(p_umap, tooltip = "text")

Интерпретация

В целом, результаты PCA и UMAP похожи. Однако UMAP сблизил некоторые страны, сделав примерно 4 кластера: 1) Наиболее богатые страны Америки (США, Канада, но также Куба и Барбадос), Океании (Австралия, Новая Зеландия), Восточной Европы (страны Балтии, РФ, Сербия, Польша), вся Западная Европа, кроме Лихтенштейна - справа внизу 2) Практически все азиатские страны, средняя азия, арабский мир, южная азия - кроме наиболее бедных; страны Кавказа, страны северной африки: как Турция, так и Морокко, а также бол-во стран латинской америки. 3) “Маленькие” cтраны: Лихтенштейн, Косово, подконтрольные маленькие территории Франции, США, и Китая, Суринам (наименьшая страна Южной Америки), Аруба (небольшой остров), маленькая Гаяна (соседка Суринама). Единственная страна, кто выбивается - это Индия, что, видимо, говорит о достаточно плохих показателях, т.к. эта страна с огромным населением и территорией, но показателями, близкими к маленьким, островным, подконтрольным другим странам государствам. Эта группа чуть слева сверху. 4) Все остальные страны Африки, а также беднейшние страны Азии и Океании, и Гаити. Группа в левом верхнем углу.

Мне кажется, UMAP сформировал крайне логичные и легко интерпретируемые кластеры, лучше, чем PCA, однако расходимость с PCA в, например, группировке африканских стран, говорит о том, что эти выделенные кластера не являются очень робастными.

11. Давайте самостоятельно увидим, что снижение размерности – это группа методов, славящаяся своей неустойчивостью. Удалите 5 случайных колонок. Проведите PCA анализ. Повторите результат 3 раза. Наблюдаете ли вы изменения в куммулятивном проценте объяснённой вариации? В итоговом представлении

данных на биплотах? С чем связаны изменения между тремя PCA?

perform_pca <- function(data) {
  pca_result <- prcomp(data, center = TRUE, scale. = TRUE)
  print(summary(pca_result))
  df_pca <- data.frame(PC1 = pca_result$x[,1], PC2 = pca_result$x[,2])
  ggplotly(ggplot(df_pca, aes(x = PC1, y = PC2, color = df$continent, text = df$Country)) + geom_point() +theme_classic(), tooltip = "text")
}

# эксперимент три раза
plots <- list()
for (i in 1:3) {
  cols_to_remove <- sample(colnames(df_scaled), 5)
  df_reduced <- df_scaled[, !(colnames(df_scaled) %in% cols_to_remove)]
  plots[[i]] <- perform_pca(df_reduced)
}
## Importance of components:
##                          PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.337 1.2586 1.15294 1.04460 0.97511 0.91890 0.85256
## Proportion of Variance 0.390 0.1131 0.09495 0.07794 0.06792 0.06031 0.05192
## Cumulative Proportion  0.390 0.5032 0.59811 0.67605 0.74397 0.80428 0.85620
##                            PC8     PC9    PC10    PC11    PC12    PC13   PC14
## Standard deviation     0.74572 0.63227 0.59430 0.54117 0.42256 0.34947 0.3325
## Proportion of Variance 0.03972 0.02855 0.02523 0.02092 0.01275 0.00872 0.0079
## Cumulative Proportion  0.89592 0.92448 0.94971 0.97063 0.98338 0.99210 1.0000
## Importance of components:
##                           PC1    PC2     PC3     PC4    PC5     PC6     PC7
## Standard deviation     2.4886 1.3437 1.14795 1.03503 0.9525 0.90151 0.70711
## Proportion of Variance 0.4424 0.1290 0.09413 0.07652 0.0648 0.05805 0.03571
## Cumulative Proportion  0.4424 0.5713 0.66546 0.74198 0.8068 0.86483 0.90054
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.63768 0.60543 0.47402 0.39184 0.35061 0.27683 0.20348
## Proportion of Variance 0.02905 0.02618 0.01605 0.01097 0.00878 0.00547 0.00296
## Cumulative Proportion  0.92959 0.95577 0.97182 0.98279 0.99157 0.99704 1.00000
## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.4471 1.4349 1.12017 1.06117 0.94571 0.87595 0.79657
## Proportion of Variance 0.4277 0.1471 0.08963 0.08043 0.06388 0.05481 0.04532
## Cumulative Proportion  0.4277 0.5748 0.66443 0.74486 0.80874 0.86355 0.90887
##                            PC8     PC9    PC10    PC11    PC12    PC13
## Standard deviation     0.66250 0.56770 0.47125 0.38386 0.32103 0.20514
## Proportion of Variance 0.03135 0.02302 0.01586 0.01052 0.00736 0.00301
## Cumulative Proportion  0.94022 0.96325 0.97911 0.98963 0.99699 1.00000
##                             PC14
## Standard deviation     2.486e-16
## Proportion of Variance 0.000e+00
## Cumulative Proportion  1.000e+00
subplot(plots)

Интерпретация

В целом, видим что общая тенденция, н-р концентрация Африканских и Европейских стран в разных частях биплотов сохраняется, однако ряд других выделенных тенденций на “полном” PCA, например очевидно разделенные кластера Африканских стран или кластер “маленьких” стран выделяются далеко не всегда. Думаю, такой метод “бутстраппинга” PCA может помочь разделить тенденции реально существующие, и появляющиеся с ходе “оверфиттинга” аналитика, смотрящего на плот. Во всех итерациях самой информативной была первая компонента, в районе 35-40% объясненной дисперсии, а 75% достигалась на 5-й компоненте.

12. Давайте самостоятельно увидим, что снижение размерности – это группа методов, славящаяся своей неустойчивостью. Создайте две дамми-колонки о том: (1) принадлежит ли страна к африканскому континенту, (2) Океании. Проведите PCA вместе с ними, постройте биплоты. Проинтерпрейтируйте результат. Объясните, почему добавление дамми-колонок не совсем корректно в случае PCA нашего типа.

df_scaled_with_dummies <- cbind(df_scaled, Is_Africa = as.numeric(df$continent == "Africa"), Is_Oceania = as.numeric(df$continent == "Oceania"))
pca_result_with_dummies <- prcomp(df_scaled_with_dummies, center = TRUE, scale. = TRUE)
loadings <- pca_result_with_dummies$rotation[, 1:2]
loadings <- data.frame(Variable = rownames(loadings), loadings)

biplot_data <- data.frame(PC1 = pca_result_with_dummies$x[,1], 
                          PC2 = pca_result_with_dummies$x[,2], 
                          Continent = df$continent,
                          Country = df$Country)

p <- ggplot() +
  geom_point(data = biplot_data, aes(x = PC1, y = PC2, color = Continent, text = Country)) +
  geom_segment(data = loadings, aes(x = 0, y = 0, xend = PC1 * scaling_factor, yend = PC2 * scaling_factor), arrow = arrow(type = "closed", length = unit(0.2, "inches")), color = "black") +
  geom_text(data = loadings, aes(x = PC1 * scaling_factor, y = PC2 * scaling_factor, label = Variable), hjust = 1.1, vjust = 1.1) +
  theme_classic()

ggplotly(p, tooltip = "text")
summary(pca_result_with_dummies)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6    PC7
## Standard deviation     2.8400 1.4851 1.3960 1.28878 1.11384 1.01306 0.9514
## Proportion of Variance 0.3841 0.1050 0.0928 0.07909 0.05908 0.04887 0.0431
## Cumulative Proportion  0.3841 0.4891 0.5819 0.66098 0.72006 0.76893 0.8120
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.92891 0.84328 0.70113 0.67720 0.62108 0.55110 0.46766
## Proportion of Variance 0.04109 0.03386 0.02341 0.02184 0.01837 0.01446 0.01041
## Cumulative Proportion  0.85312 0.88699 0.91040 0.93223 0.95060 0.96506 0.97548
##                           PC15    PC16    PC17    PC18    PC19    PC20
## Standard deviation     0.39263 0.35858 0.34148 0.26517 0.20133 0.06900
## Proportion of Variance 0.00734 0.00612 0.00555 0.00335 0.00193 0.00023
## Cumulative Proportion  0.98282 0.98894 0.99449 0.99784 0.99977 1.00000
##                             PC21
## Standard deviation     2.959e-16
## Proportion of Variance 0.000e+00
## Cumulative Proportion  1.000e+00

Интерпретация

По идее, PCA предполагает линейные взаимосвязи между переменными и лучше всего работает с непрерывными переменными. Поэтому добавление бинарных данных теоретически некорректно, так как PCA будет пытаться интерпретировать эти бинарные значения как непрерывные. Однако у меня, несмотря на то, что соответствующие собственные вектора отображаются на плоте и саммари результатов явно обновилось, положение точек практически не изменилось на биплоте. Видимо, вклад других переменных в первые две компоненты крайне велик так, что добавление бинарных мало влияет. Возможно, если бы мы добавили категориальную переменную, которая вносила бы значительный вклад в первые компоненты, это бы привело к искажению данных.